Revisiting the slow dynamics of a silica melt using Monte Carlo simulations 
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We implement a standard Monte Carlo algorithm to study the slow, equilibrium dynamics of a 
silica melt in a wide temperature regime, from 6100 K down to 2750 K. We find that the average 
dynamical behaviour of the system is in quantitative agreement with results obtained from molec- 
ular dynamics simulations, at least in the long-time regime corresponding to the alpha-relaxation. 
By contrast, the strong thermal vibrations related to the Boson peak present at short times in 
molecular dynamics are efficiently suppressed by the Monte Carlo algorithm. This allows us to 
reconsider silica dynamics in the context of mode-coupling theory, because several shortcomings 
of the theory were previously attributed to thermal vibrations. A mode-coupling theory analysis 
of our data is qualitatively correct, but quantitative tests of the theory fail, raising doubts about 
the very existence of an avoided singularity in this system. We discuss the emergence of dynamic 
heterogeneity and report detailed measurements of a decoupling between translational diffusion and 
structural relaxation, and of a growing four-point dynamic susceptibility. Dynamic heterogeneity 
appears to be less pronounced than in more fragile glass-forming models, but not of a qualitatively 
different nature. 

PACS numbers: 02.70.Ns, 64.70.Pf, 05.20.Jj 



I. INTRODUCTION 

Numerical simulations play a major role among stud- 
ies of the glass transition since, unlike in experimental 
works, the individual motion of a large number of parti- 
cles can be followed at all times Computer simula- 
tions usually study Newtonian dynamics (ND) by solving 
a discretized version of Newton's equations for a given in- 
teraction between particles . Although this is the most 
appropriate dynamics to study molecular liquids, it can 
be interesting to consider alternative dynamics that are 
not deterministic, or which do not conserve the energy. In 
colloidal glasses and physical gels, for instance, the par- 
ticles undergo Brownian motion arising from collisions 
with molecules in the solvent, and a stochastic dynamics 
is more appropriate Q . Theoretical considerations might 
also suggest the study of different sorts of dynamics for a 
given interaction between particles, for instance to assess 
the role of conservation laws [1, Q, H] and structural in- 
formation [1, 01- Of course, if a given dynamics satisfies 
detailed balance with respect to the Boltzmann distri- 
bution, all structural quantities remain unchanged, but 
the resulting dynamical behaviour might be very differ- 
ent. In this paper, we study the relaxation dynamics of a 
commonly used theoretical model for silica using Monte 
Carlo simulations and compare the results with previous 
ND studies for both the averaged dynamical behaviour 
and the spatially heterogeneous dynamics of this system. 

Several papers have studied in detail the influence of 
the chosen microscopic dynamics on the dynamical be- 
haviour in a simple glass-former, namely a binary mix- 
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ture of Lennard-Jones particles [H, Q. Gleim, Kob and 
Binder Q studied Stochastic Dynamics where a friction 
term and a random noise are added to Newton's equa- 
tions, the amplitude of both terms being related by a 
fluctuation-dissipation theorem. Szamel and Flenner [lO| 
used Brownian dynamics, in which there are no momenta, 
and positions evolve with a Langevin dynamics. Berthier 
and Kob [ll| employed Monte Carlo dynamics, where the 
potential energy between two configurations is used to 
accept or reject a trial move. The equivalence between 
these three stochastic dynamics and the originally stud- 
ied ND was established at the level of the averaged dy- 
namical behaviour, except at very short times where dif- 
ferences are indeed expected. However, important differ- 
ences were found when dynamic fluctuations were consid- 
ered, even in the long-time regime comprising the alpha- 
relaxation [i,l3,IIH- 

Silica, the material studied in the present work, is dif- 
ferent from the previously considered Lennard-Jones case 
in many aspects which all motivate our Monte Carlo 
study of the Beest, Kramer and van Santen (BKS) model 
for silica [l^. First, the BKS model was devised to rep- 
resent a real material, making our conclusions more di- 
rectly applicable to experiments. Second, the temper- 
ature evolution of relaxation times is well-described by 
a simple Arrhenius law at low temperatures, typical of 
strong glass-formers, which are commonly believed to be- 
long to a somewhat different class of materials, so that 
qualitative differences might be expected with more frag- 
ile, super- Arrheniusly relaxing materials. Third, the on- 
set of slow dynamics in fragile materials is often said to be 
accurately described by the application of mode-coupling 
theory, at least over an intermediate window of 2 to 3 
decades of relaxation times Mode-coupling theory 
(MCT) formulates in particular a series of quantitative 
predictions regarding the time, spatial, and temperature 
dependences of dynamic correlators. In the case of sil- 
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ica melts, previous analysis reported evidence in favour 
of a narrower temperature regime where MCT can be 
applied, but the test of several theoretical predictions 
was either seriously affected, or even made impossible 
by the presence of strong short-time thermal vibrations 
related to the Boson peak in this material [T3|- These 
vibrations affect the time dependence of the correlators 
much more strongly in silica than in Lcnnard- Jones sys- 
tems, which constitutes a fourth difference between the 
two systems. Using Monte Carlo simulations we shall 
therefore be able to revisit the MCT analysis performed 
in Refs. [l^. Fifth, while detailed analysis of dynamic 
heterogeneity is available for fragile materials, a com- 
paratively smaller amount of data is available for strong 
materials 0, 0, EB] , and we shall therefore investigate 
issues that have not been addressed in previous work. 

The paper is organized as follows. In Sec. [IT] we give 
details about the simulation technique and compare its 
efficiency to previously studied dynamics. In Sec. IIIII we 
present our numerical results about the averaged dynam- 
ics of silica in Monte Carlo simulations, while Sec. IIVI 
deals with aspects related to dynamic heterogeneity. Fi- 
nally, Sec. |V] concludes the paper. 

II. SIMULATING SILICA USING MONTE 
CARLO DYNAMICS 

Our aim is to study a non-Newtonian dynamics of the 
glass- former silica, Si02. Therefore, we must first choose 
a reliable model to describe the interactions in this two- 
component system made of Si and O atoms, and then 
design a specific stochastic dynamics which we require 
to be efficient and to yield the same static properties as 
Newtonian dynamics. 

Various simulations have shown that a reliable pair 
potential to study silica in cornputer simulations is the 
one proposed by BKS [H, [H, Q M, [H, [13, [3- The 
functional form of the BKS potential is 

W = — ^ + A^pexpi-Baisr) - (1) 

where a, /? € [Si, O] and r is the distance between the 
atoms of type a and /3. The values of the constants 
Qa, qp, ^ap, Bai3, and Caf3 cau be found in Ref. For 
the sake of computational efficiency the short range part 
of the potential was truncated and shifted at 5.5 A. This 
truncation also has the benefit of improving the agree- 
ment between simulation and experiment with regard to 
the density of the amorphous glass at low temperatures. 
The system investigated has A^gi = 336 and A^o = 672 
atoms in a cubic box with fixed size L = 24.23 A, so 
that the density is p = 2.37 g/cm^, close to the exper- 
imentally measured density at atmospheric pressure of 
2.2 g/cm3 [3. 

Once the pair interaction is chosen, we have to de- 
cide what stochastic dynamics to implement. Previ- 
ous studies in a Lennard- Jones system concluded that 



among Monte Carlo (MC), Stochastic Dynamics (SD) 
and Brownian Dynamics (BD), MC was by far the most 
efficient algorithm because relatively larger incremen- 
tal steps can be used while maintaining detailed bal- 
ance, which is impossible for SD and BD where very 
small discrctized timcstcps are needed to maintain the 
fluctuation-dissipation relation between noise and fric- 
tion terms [Tlj. Given the generality of this argument, it 
should carry over to silica, and we decided to implement 
MC dynamics to the BKS model. An additional justifi- 
cation for our choice stems from unpublished results by 
Horbach and Kob who performed preliminary investiga- 
tions of the SD of BKS silica [2^ . Using a friction term 
similar in magnitude to the one used in Lennard- Jones 
simulations was however not enough to efficiently sup- 
press short-time clastic vibrations. Using an even larger 
friction term would probably damp these vibrations, but 
would also make the simulation impractically slow. 

A standard Monte Carlo dynamics for the pair po- 
tential in Eq. ([1]) should proceed as follows. In an ele- 
mentary MC move, a particle, i, located at the position 
Ti is chosen at random. The energy cost, AEi, to move 
particle i from position to a new, trial position -I- 6r 
is evaluated, Sr being a random vector comprised in a 
cube of linear length <5niax centered around the origin. 
The Metropolis acceptance rate, p = min(l, e^0^Ei/kB -j^ 
where P = 1/T is the inverse temperature, is then 
used to decide whether the move is accepted or rejected. 
In the following, one Monte Carlo timestep represents 
N = Nsi + No attempts to make such an elementary 
move, and timescales are reported in this unit. Temper- 
ature will be expressed in Kelvin. Monte Carlo simula- 
tions can of course be made even more efficient by imple- 
menting for instance swaps between particles, or using 
parallel tempering. The dynamical behaviour, however, 
is then strongly affected by such non-physical moves and 
only equilibrium thermodynamics can be studied. Since 
we want to conserve a physically realistic dynamics, we 
cannot use such improved schemes. 

An additional difficulty with Eq. ([1]) as compared to 
Lennard-Jones systems is the Coulombic interaction in 
the first term. Such a long-range interaction means that 
the evaluation of the energy difference AEi needed to 
move particle i in an elementary MC step requires a sum 
over every particles j ^ i, and over their repeated im- 
ages due to periodic boundary conditions. Of course 
the sum can be efficiently evaluated using Ewald summa- 
tions techniques, as is commonly employed in ND simu- 
lations 0. We note, however, that Ewald techniques are 
better suited for ND than for MC since in ND the po- 
sitions and velocities of all particles are simultaneously 
updated so that the Ewald summation is performed once 
to update all particles. In MC simulations, each single 
move requires its own Ewald summation, and this re- 
mains computationally very costly. 

For the BKS potential in Eq. ([1]) it was recently shown 
that a simple truncation can be performed which makes 
the range of the Coulombic interaction term finite [2l| . 
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FIG. 1: Self-intermediate scattering function for silicon, 
Eq. (O, at T = 4000 K and |q| = 1.7 A"^ for various values 
of (Smax. Inset: The evolution of the relaxation time with 5max 
unambiguously defines an optimal value (5max ~ 0.65 A for 
efficient Monte Carlo simulations. 



Detailed ND simulations have shown that in the range 
of temperatures presently accessible to computer exper- 
iments, no difference can be detected between the finite 
range and the infinite range versions of the BKS poten- 
tial for a wide variety of static and dynamic properties. 
Therefore we build on this work and make the replace- 
ment m 
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c), for r < Tc 



(2) 



while 1/r — > for r > Tc. This amounts to smoothly 
truncating the potential at a finite range, Tc, maintaining 
both energy and forces continuous at the cutoff r = Tc- 
The physical motivation for this form of the trunca- 
tion was given by Wol f I23 |. and discussed in several 
more recent papers [2l|, |23||. Following Ref . [2l| , we fix 
Tc = 10.14 A. Once the potential is truncated, MC sim- 
ulations become much more efficient, and much simpler 
to implement. Furthermore, this will allow us to perform 
detailed comparisons of the dynamics of BKS model of 
silica where the Wolf truncation is used for both ND and 
MC in the very same manner, so that any difference be- 
tween the two sets of data can be safely assigned to the 
change of microscopic dynamics alone, while reference to 
earlier work done using Ewald summations is still quan- 
titatively meaningful. 

The one degree of freedom that remains to be fixed 
is f^max, which determines the average lengthscale of el- 
ementary moves. If chosen too small, energy costs are 
very small and most of the moves are accepted, but the 
dynamics is very slow because it takes a long time for par- 
ticles to diffuse over the long distances needed to relax 
the system. On the other hand too large displacements 
will on average be very costly in energy and acceptance 
rates can become prohibitively small. We seek a com- 
promise between these two extremes by monitoring the 
dynamics at a moderately low temperature, T = 4000 K, 



for several values of (5max • As a most sensitive indicator of 
the relaxational behaviour, we measure the contribution 
from the specie a [a = Si, O) to the self-intermediate 
scattering function defined by 
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We make use of rotational invariance to spherically av- 
erage over wave vectors of comparable magnitude, and 
present results for |q| = 1.7 A"-'^, which is the location of 
the pre-peak observed in the static structure factor S{q) 
of the liquid. This corresponds to the typical (inverse) 
size of the Si04 tetrahedra. In Fig. [T] we present our re- 
sults for ^max values between 0.3 and 1.0 A. As expected 
we find that relaxation is slow both at small and large 
values of (5niax, and most efficient for intermediate values. 
Interestingly we also note that the overall shape of the 
self-intermediate scattering function does not sensitively 
depend on (5niax over this wide range. We can therefore 
safely fix the value of (5max based on an efficiency cri- 
terium alone. 

We define a typical relaxation time, Tq,, as 



(4) 



and show its (Jmax dependence in the inset of Fig. [T] 
A clear minimum is observed at the optimal value of 
(^max ~ 0.65 A, which we therefore use throughout the 
rest of this paper. This distance corresponds to a squared 
displacement of 0.4225 A^, which is very close to the 
plateau observed at intermediate times in the mean- 
squared displacements (see Fig. [2] below). This plateau 
can be taken as a rough measurement of the "cage" size 
for the particles, so that MC simulations are most effi- 
cient when the cage is most quickly explored. This ar- 
gument and the data in Fig. [2] suggest that the location 
of the minimum should only be a weak function of tem- 
perature, but we have not verified this point in detail. 
Therefore we keep the value of (5max constant at all tem- 
peratures. An alternative would be to optimize it at each 
T and then carefully rescale timescales between runs at 
different temperatures. 

What about the relative efficiency between MC and 
ND? If we compare the relaxation measured at T = 
4000 K, we find ~ 13400 Monte Carlo steps, while 
Ta ^ 4.7 ps for ND. When using a discretized timestep of 
1.6 fs, this means that, when counting in number of inte- 
gration timestcps, MC dynamics is « 5 times slower than 
ND. This result contrasts with the results obtained in a 
Lennard-Jones mixture where MC dynamics was about 
2 times faster than ND [ll|. We attribute this relative 
loss of efficiency to the existence of strong bonds between 
Si and O atoms in silica, which have no counterparts in 
Lennard-Jones systems. It is obvious that strong bonds 
are very hard to relax when using sequential Monte Carlo 
moves, as recently discussed in Ref. 24 1. 
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We have performed simulations at temperatures be- 
tween T = 6100 K and T = 2750 K, the latter be- 
ing smaller than the fitted mode-coupling temperature, 
Tc = 3330 K For each temperature we have sim- 

ulated 3 independent samples to improve the statis- 
tics. Initial configurations were taken as the final con- 
figurations obtained from previous work performed with 
ND [U , so that production runs could be started imme- 
diately. For each sample, production runs lasted at least 
IStq (at T ~ 2750 K), much longer for higher tempera- 
tures, so that statistical errors in our measurements are 
fairly small. We have performed a few runs for a larger 
number of particles, namely N = 8016 particles, to inves- 
tigate finite size effects which are known to be relevant 
in silica [l^, [H, H^l , and the results will be discussed in 

secnni 



III. ANALYSIS OF AVERAGED TWO-TIME 
CORRELATORS 

In this section we report our results about the time be- 
haviour of averaged two-time correlators, we compare the 
Monte Carlo results to Newtonian dynamics, and we per- 
form a quantitative mode-coupling analysis of the data. 



A. Intermediate scattering function and 
mean-squared displacements 

The self- intermediate scattering function, Eq. ([3|), is 
shown in Fig. [2] for temperatures decreasing from T = 
6100 K down to T = 2750 K for Si and O atoms at 
Iql = 1.7 A-\ These curves present well-known features. 
Dynamics at high temperature is fast and has an expo- 
nential nature. When temperature is decreased below 
T 4500 K, a two-step decay, the slower being strongly 
non-cxponcntial, becomes apparent. Upon decreasing 
the temperature further, the slow process dramatically 
slows down by about 4 decades, while clearly conserv- 
ing an almost temperature-independent, non-exponential 
shape, as already reported for ND [l^ . 

We also find that the first process, the decay towards a 
plateau, slows down considerably when decreasing tem- 
perature, although less dramatically than the slower pro- 
cess. The fastest process, called 'critical decay' in the 
language of mode-coupling theory 0], is not observed 
when using ND, because it is obscured by strong ther- 
mal vibrations occurring at high frequencies (in the THz 
range). Clearly, no such vibrations are detected in the 
present results which demonstrates our first result: MC 
simulations very efficiently suppress the high-frequency 
oscillations observed with ND. 

Although the plateau seen in Fs{q,t) is commonly in- 
terpreted as 'vibrations of a particle within a cage', the 
data in Fig. [2] discard this view. From direct visualization 
of the particles' individual dynamics it is obvious that vi- 
brations take place in just a few MC timestcps, while the 
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FIG. 2: Top: Self-intermediate scattering function, Eq. ((2]), 
for lq| = 1.7 and temperatures T = 6100, 4700, 4000, 
3580, 3200, 3000, and 2750 K (from left to riglit). Bottom: 
Mean-squared displacement, Eq. ([5]), for the same tempera- 
tures in the same order. 



decay towards the plateau can be as long as 10"* time 
units at the lowest temperatures studied here. This de- 
cay is therefore necessarily more complex, most probably 
cooperative in nature. This interpretation is supported 
by recent theoretical studies where a plateau is observed 
in two-time correlators of lattice models where local vi- 
brations arc indeed completely absent [1^. A detailed 
atomistic description of this process has not yet been re- 
ported, but would indeed be very interesting. 

Next, wc study the mean-squared displacement defined 
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as 



^(|r,(i)-r,;(0)|^), 



(5) 



and wc present its temperature evolution in Fig. [21 for 
both Si and O atoms. The evolution of A^r(t) mirrors 
that of the self-intermediate scattering function, and the 
development of a two-step relaxation process is clear from 
these figures. Because we are studying a stochastic dy- 
namics, displacements are diffusive at both short and 
long timescales. This constitutes an obvious, expected 
difference between ND and MC simulations: data clearly 
cannot match at very small times. The goal of the present 
study is therefore to determine whether the dynamics 
quantitatively match at times where the relaxation is not 
obviously ruled by short-time ballistic/diffusive displace- 
ments. 

The plateau observed in Fs{q,t) now translates into 
a strongly sub-diffusive regime in the mean-squared dis- 
placements separating the two diffusive regimes. At the 
lowest temperature studied, when t changes by three 
decades from 2 x 10^ to 2 x 10^, the mean-squared dis- 
placement of Si changes by a mere factor 4.6 from 0.16 to 
0.074. Particles are therefore nearly arrested for several 
decades of times, before eventually entering the diffusing 
regime where the relaxation of the structure of the liquid 
takes place. 



B. Comparison to Newtonian dynamics 

The previous subsection has shown that the Monte 
Carlo dynamics of silica is qualitatively similar to the one 
reported for ND, apart at relatively short times where 
the effect of thermal vibrations is efficiently suppressed 
and the dynamics is diffusive instead of ballistic. We 
now compare our results more quantitatively with the 
dynamical behaviour observed using ND. 

To this end, we compare first the temperature evolu- 
tion of the relaxation times, Ta{T), defined in Eq. in 
Fig. [31 Here, we use a standard representation where an 
Arrhenius slowing down over a constant energy barrier 
E, with an attempt frequency I/tq, 



ToCXp 



E 



(6) 



appears as a straight line. To compare both sets of data 
we rescale the ND data by a common factor, = 0.31 fs, 
which takes into account the discretization timestep and 
the efficiency difference discussed in the previous section; 
to will be kept constant throughout this paper. We find 
that the temperature evolution of the alpha-relaxation 
time measured in MC simulations is in complete quan- 
titative agreement with the one obtained from ND, over 
the complete temperature range. This proves that Monte 
Carlo techniques can be applied not only to study static 
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FIG. 3: Temperature evolution of the alpha-relaxation time 
Ta{T) for silicon (squares) and oxygen (circles), and inverse 
self-diffusion constant for silicon (up triangles) and oxygen 
(down triangles), vertically shifted for clarity. Open symbols 
are for ND (times rescaled by to ~ 0.31 fs) closed symbols 
for MC. Full lines are Arrhenius fits below T « 3700 K with 
activation 5.86, 5.60, 5.43, and 4.91 eV (from top to bottom). 
An Arrhenius fit for high temperatures is also presented for 
D(0) with activation energy 2.76 eV. The dashed line is a 
power law fit, ~ (T - Tc)~^ , with Tc = 3330 K and 7 = 
2.35. 



properties of silica, but also its long-time dynamic prop- 
erties. 

In Fig. [3] we also show the temperature evolution of the 
self-diffusion constant, defined from the long-time limit 
of the mean-squared displacement as 



D = lim 



6t 



(7) 



The behaviour of the (inverse) diffusion constant is quali- 
tatively very close to the one of the alpha-relaxation time, 
and we again find that ND and MC dynamics yield re- 
sults in full quantitative agreement. 

As expected for silica, we find that at low tempera- 
tures below T w 3700 K, relaxation timescales and diffu- 
sion constant change in an Arrhenius fashion described 
by Eq. ([6]). We find, however, that the observed acti- 
vation energies display small variations between differ- 
ent observables, from 5.86 eV for Ta(Si) to 4.91 eV for 
1 / D(0) ■ These values compare well with previous analy- 
sis [l3| , and with experimental findings [29| . 

From Fig. [31 it is clear that Arrhenius behaviour is 
obeyed below T « 3700 K only, while the data bend up 
in this representation for higher temperatures. This be- 
haviour was interpreted in terms of a fragile to strong 
behaviour of the relaxation timescales in several pa- 
pers [H, [10, Hll, despite the fact that fragility is usu- 
ally defined experimentally by considering data on a 
much wider temperature window close to the experimen- 
tal glass transition. To rationalize these findings, Hor- 
bach and Kob analyzed the data using mode-coupling 
theory predictions [l3|. In particular they suggest to fit 
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FIG. 4: Decoupling data for oxygen and silicon. We plot the 
product Dra taken from the data shown in Fig. |3] and nor- 
malize the product by its value at To = 4700 K such that 
deviations from 1 indicates non-zero decoupling. Open sym- 
bols are for ND, closed symbols for MC. Decoupling is similar 
for both types of dynamics. 



the temperature dependence of Tq as 



(8) 



with Tc « 3330 K and 7 « 2.35. This power law fit is also 
presented in Fig. [3] as a dashed line. Its domain of validity 
is of about 1 decade, which is significantly less than for 
more fragile materials with super- Arrhenius behaviour of 
relaxation timescales Q . 

It is interesting to note that a simpler interpretation of 
this phenomenon could be that this behaviour is nothing 
but a smooth crossover from a non-glassy, homogeneous, 
high-temperature behaviour to a glassy, heterogeneous, 
low temperature behaviour, as found in simple models of 
strong glass- forming liquids [3^ . In Fig. [31 we implement 
this simpler scenario by fitting high temperature data 
with an Arrhenius law, as is sometimes done in the anal- 
ysis of experimental data (33l | . Such a fit works nicely for 
high temperatures, from T = 6100 to 4700 K, but breaks 
down below T « 4000 K. A physical interpretation for 
this high-temperature Arrhenius behaviour was offered 
in Ref. [sj. This shows that analyzing silica dynamics 
in terms of a simple crossover occurring around 4000 K 
between two simple Arrhenius law is indeed a fair descrip- 
tion of the data which does not require invoking a more 
complex fragile to strong crossover being rationalized by 
the existence of an avoided mode-coupling singularity. 

The difference found above for the activation energies 
describing Tq, and 1 /D for both species implies that these 
quantities, although both devised to capture the tem- 
perature evolution of single particle displacements have 
slightly different temperature evolutions and are not pro- 
portional to one another. This well-known feature im- 
plies the existence of a "decoupling" between transla- 
tional diffusion and structural relaxation in silica, as doc- 
umented in previous papers [13. In Fig. a we report 
the temperature evolution of the product D{T)Ta{q,T) 
which is a pure constant for a simply diffusive particle 



where Ta(q,T) = l/{q^D). We normalize this quantity 
by its value at To — 4700 K, so that any deviations from 1 
indicates a non-zero decoupling (35l . ISq . As expected we 
find that the product is not a constant, but grows when 
temperature decreases. Remarkably, although this quan- 
tity is a much more sensitive probe of the dynamics of the 
liquid, its temperature evolution remains quantitatively 
similar for both ND and MC dynamics. This shows that 
equivalence of the dynamics between the two algorithms 
holds at the level of the complete distribution of particle 
displacements, even for those tails that are believed to 
dictate the observed decoupling. 

In Sec. llVi we shall explore in more detail the het- 
erogeneous character of the dynamics of silica, closely 
related to the decoupling discussed here. It is however 
interesting to try and infer the amount of decoupling pre- 
dicted for silica at temperatures close to the experimen- 
tal glass transition, Tg w 1450 K. The glass transition 
temperature of BKS silica deduced from extrapolation of 
viscosity measurements is close to the experimental one, 
T^KS ^ 1350 K [H. Extrapolating the data in Fig. g] 
down to 1400 K predicts a decoupling of about 40 for 
Si dynamics, about 7 for O dynamics. The difference 
between Si and O dynamics was recently explained in 
Ref. [131, where it was noted that oxygen diffusion is in 
fact possible with no rearrangement of the tetrahedral 
structure of silica involved. Moreover, it is interesting to 
note that the amount of decoupling found here is smaller 
than experimental findings in fragile materials close to 
their glass transition [37[, but is nonetheless clearly dif- 
ferent from zero. This suggests that even strong materi- 
als display dynamically heterogeneous dynamics, but its 
effect seems less pronounced than in more fragile mate- 
rials. 

Theoretically, an identical temperature evolution of 
the alpha relaxation timescale for MC and ND is an im- 
portant prediction of mode-coupling theory [3] because 
the theory uniquely predicts the dynamical behaviour 
from static density fluctuations. Gleim et al. argue that 
their finding of a quantitative agreement between SD and 
ND in a Lennard- Jones mixture is a nice confirmation 
of this non-trivial mode-coupling prediction Q. Szamel 
and Flenner [lo| confirmed this claim using BD, and ar- 
gued further that even deviations from mode-coupling 
predictions are identical, a statement that was extended 
to below the mode-coupling temperature by Berthier and 
Kob [ll[ . In the present work we extend these findings to 
the case of silica over a large range of temperatures, which 
goes far beyond the temperature regime where MCT can 
be applied. Therefore, we conclude that such an inde- 
pendence of the glassy dynamics of supercooled liquids to 
their microscopic dynamics, although predicted by MCT, 
certainly has a much wider domain of validity than the 
theory itself. Finally, we note that the deviations from 
MCT predictions observed in Fig. [3] cannot be attributed 
to coupling to currents which are expressed in terms of 
particle velocities. In our MC simulations we have no 
velocities, so that avoiding the mode-coupling singular- 
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FIG. 5: Self-intermediate scattering function for fixed T and 
q — 1.7 A^^, obtained in MC and MD simulations for two 
system sizes. The time axis in MD data is rescaled by to = 
0.31 fs to obtain maximum overlap with MC results, and the 
same factor is used for the two sizes. Larger systems relax 
faster and the amplitude of this finite size effect is the same 
for both dynamics. 



ity is not due to the hydrodynamic effects pointed out 
in Ref. Q (see Ref. [1^ for more recent theoretical view- 
points) . 

The last comparison to ND wc want to discuss con- 
cerns the study of finite size effects. It was shown that 
the long-time dynamics of silica is fairly sensitive to sys- 
tem size, and there are detectable differences when the 
number of particles is changed from 1000 to 8000 [H,[23- 
Such a large effect is not observed in more fragile ma- 
terials [1]. It was suggested that short-time thermal 
vibrations, stronger in silica than in simpler models, 
are responsible for this system size dependence [13, . 
Therefore, it could be expected that by efficiently sup- 
pressing these vibrations finite size effects should be re- 
duced. But this is not what happens. In Fig. O we 
show self-intermediate scattering functions measured at 
T = 3580 K and |q| = 1.7 A"! in both ND and MC 
for two system sizes, = 1008 and = 8016 par ticles. 
Such data have been presented for ND before [23|, and 
our results agree with these earlier data. The amplitude 
of the vibrations observed for t/to = ^ ~ 10"^ is smaller 
and the long time dynamics is faster when N is larger. 
For MC we find that high-frequency vibrations and the 
corresponding finite size effects are indeed suppressed, 
but the finite size effect for long-time relaxation, some- 
what surprisingly, survives in our MC simulations, and 
can therefore not be attributed to high-frequency thermal 
vibrations. Recent studies of the vibration spectrum and 
elastic properties at T = of amorphous media have sug- 
gested the existence of large-scale structures [1^ : these 
objects are potential candidates to account for the size 
effect found at long times. It should then be explained 
how these spatial structures affect the long time dynam- 
ics, and why a finite size simulation box at the same 
time affects the absolute value of the alpha-relaxation 
timescale but leaves unchanged many of its detailed prop- 
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C. Mode-coupling analysis of dynamic correlators 

We now turn to a more detailed analysis of the 
shape and wave vector dependences of two-time correla- 
tion functions, revisiting in particular the mode-coupling 
analysis performed by Horbach and Kob in Refs. p^ . 
They argue that MCT can generally be applied to de- 
scribe their silica data, and attribute most of the devia- 
tions that they observe to short-time thermal vibrations 
supposedly obscuring the "true" MCT behaviour. We 
are therefore in a position to verify if their hypothesis is 
correct. 

When applied to supercooled liquids, MCT formulates 
a scries of detailed quantitative predictions regarding the 
time, wave vector, and temperature dependences of two- 
time dynamical correlators close to the mode-coupling 
singularity. In particular, MCT predicts that correlation 
functions should indeed decay in the two-step manner re- 
ported in Fig. [21 Moreover, for intermediate times corre- 
sponding to the plateau observed in correlation functions, 
an approximate equation can be derived which describes 
the correlator close enough to the plateau Q. The fol- 
lowing behaviour is then predicted. 



Fs{q,t)^ fg + hgF{t), 



(9) 



where F{t) is the so-called /3-correlator which is inde- 
pendent of the wave vector, and whose shape depends on 
a few parameters: the reduced distance from the mode- 
coupling temperature, e = \T — Tc\/Tc, and a parame- 
ter describing the MCT critical exponents, A. Once A 
is known various exponents (0,6,7) ^-rc known, which 
describe, in particular, the short-time behaviour of F(t) 
when Fs{q, t) approaches the plateau, F{t) ~ t~°, and its 
long-time behaviour when leaving the plateau, F{t) ~ t^. 
The exponent 7 was introduced in Eq. ^ and describes 
the temperature evolution of the relaxation time . 

Several properties follow from Eq. ([9]). If one works at 
fixed temperature and varies the wave vector, the follow- 
ing quantity. 



0(t) - m m - F{t') 



(f){t") - (f){t') F{t")-F{t'Y 



(10) 



where (j)(t) stands for a two-time correlation function, 
should become independent of q. In Eq. ((TU]), t' and t" 
are two arbitrary times taken in the plateau regime. This 
is called the "factorization property" in the language of 
MCT. We follow Ref. [H and show in Fig. [6] the func- 
tion R{t) in Eq. (jlOp using self-intermediate scattering 
functions for different q and for different species (Si and 
O) at fixed temperatures, T = 3580 K and T = 3000 K, 
choosing times comparable to those reported in Ref. [l^ , 
namely t" = 82 and t' ^ 760 for T 3580 K, and 
t" = 1360 and t' = 66700 for T = 3000 K. Ahhough 
the factorization property seemed to hold quite well in 




t 

FIG. 6: Test of the factorization property, Eq. (|10|) using 
Fs{q,t) from Si and O dynamics, and wave vectors between 
0.8 and 4 A~\ for T = 3580 K and 3000 K. The data do not 
show collapse for times t' < t < t" , and factorization does not 
work very well. 



the ND data, this is no more the case for our MC data, 
and R{t) retains a clear q dependence between t' and t": 
no collapse of R{t) can be seen in the regime t' < t < t" 
in Fig. [51 The reason is clear from Fig. O due to ther- 
mal vibrations, the intermediate plateau was very flat in 
ND, but it has much more structure in our MC data. 
It was therefore easier to collapse the ND data in this 
regime than the present MC data for which a better 
agreement might have been expected. In the case of 
the factorization property, the presence of thermal vi- 
brations in fact favours a positive reading of the data, 
which become much less convincing when these vibra- 
tions are suppressed. Gleim and Kob had reached an 
opposite conclusion in the case of a Lennard-Joncs sys- 
tem They found that suppressing vibrations made 
the mode-coupling analysis of the beta-relaxation more 
convincing, suggesting that MCT describes the Lennard- 
Joncs system more accurately than silica. 

Next we perform a test of the theory which had not 
been possible with ND data. We investigate in detail if 
the behaviour predicted by Eq. ^ is correct for both 
short and long times. This test is not possible using 
ND because the approach to the plateau is mainly ruled 
by thermal vibrations (see for instance the ND data pre- 
sented in Fig.[n|). In Fig.[7]we show that a "critical decay" 



FIG. 7: Self-intermediate scattering function at fixed T = 
3580 K and various wave vectors, q = 0.8, 1.2, 1.7, 2.4, 3.2, 
and 4 (from right to left). Dashed lines show fits at 

intermediate times using Eq. The inset shows the q- 

dependence of the fitting parameters hq and fq. Note that 
the time domains over which the fits apply shift with q. 



does indeed show up when thermal vibrations are over- 
damped and no oscillations can be seen. To check in more 
detail if this behaviour is indeed in quantitative agree- 
ment with the MCT predictions, we fit the Fs{q,t) data 
at T = 3580 K, i.e. slightly above = 3330 K, for sev- 
eral wave vectors q using the /3-correlator obtained from 
numerical integration of the mode-coupling equation. To 
get the fits shown in Fig. [7| we have to fix the distance to 
the mode-coupling temperature e and the value A = 0.71 
both taken from Ref. [ij], and yielding a = 0.32, b = 0.62 
and 7 = 2.35. Additionally we have to adjust the micro- 
scopic timescale. Moreover, for each wave vector we have 
to fix hq and fq which respectively correspond to the am- 
plitude of the /?-correlator and the height of the plateau 
in Fs{q,t). Finally, there are two additional "hidden" 
free parameters in each of these fits: the somewhat arbi- 
trarily chosen boundaries of the time domain where the 
fitting function describes the data. We then get the fits 
shown with dashed lines in Fig. [71 which are of a quality 
comparable to the ones usually found in the MCT liter- 
ature 0. The parameters hq and fq are also shown in 
the inset of Fig. [71 and behave qualitatively as in similar 
studies. Inspection of Fig. [71 reveals that the use of such 
freedom to fit the data allows a qualitatively correct de- 
scription of the data, although clearly the time domain 
over which each wave vector is fit systematically shifts 
when q changes, and we could not simultaneously fit the 
data at both small and large q by fixing the time inter- 
val of the fit. This failure is consistent with the above 
finding that the factorization property is not satisfied. 

Therefore we conclude that MCT provides a quali- 
tatively correct description of our data in the plateau 
regime, with no satisfying quantitative agreement, even 
in the absence of short-time thermal vibrations. One has 
therefore to argue that the data are taken too far from 
the transition for MCT to quantitatively apply to silica. 
However, since it is not possible to get data closer to 
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fg{T) = f,{T,) + a^/T~T, 



T<Tr.. 



(12) 
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while fq{T > Tc) ~ fq(Tc). The non-analytic behaviour 
of fq at Tc is a further characteristic feature of the mode- 
coupling singularity. Since we can easily take data for 
T < Tc which are arguably not influenced by thermal 
vibrations, we can directly check for the presence of the 
square-root singularity, Eq. (|12p . This is done in the bot- 
tom panel of Fig. [51 where we also show data obtained 
from ND simulations. That the latter are strongly in- 
fluenced by thermal vibrations is clear, since they sys- 
tematically lie below the MC data and have a stronger 
temperature dependence close to Tc- However, even the 
MC data clearly indicate that fq{T) is better described 
by a non-singular function of temperature, compatible 
with the simple linear behaviour expected to hold at very 
low temperatures. The temperature dependence of the 
plateau height therefore explains why time temperature 
superposition does not hold in the late /3-regime, but the 
linear temperature behaviour indicates that there is no 
clear sign, from our data, of the existence of a "true" 
underlying singularity at Tc- 



FIG. 8: Top: Test of time-temperature superposition, 
Eq. mj. The dashed line is a stretched exponential func- 
tion with /3 = 0.87. Superposition holds at large rescaled 
times, but fails in the /3-regime because the plateau height in- 
creases when T decreases. Bottom: extracted plateau height 
as a function of temperature fitted with a linear dependence 
(full line) and with a square root singularity, Eq. (|12|l (dashed 
line). Open symbols are for ND, closed symbols for MC. No 
singular behaviour of fq is visible in either set of data. 



the transition (recall that the transition does not exist), 
the domain of validity of the theory then would become 
vanishingly small. 

We now turn to longer timescales and show in Fig. [5] 
a test of the time-temperature superposition prediction 
of the theory which states that correlators at fixed q but 
different temperatures should scale as 



F.,iq,t)^Tq 



Ta{q) 



(11) 



where ^q{x) ~ exp(— x'^'-'?-') and for times in the a- 
regime. When high temperatures outside the glassy 
regime are discarded Eq. (jlip works correctly when the 
scaling variable t/Ta is not too small, but fails more 
strongly in the late /3-rcgime. Scaling in the /3-regime 
is often one of the most successful prediction of MCT, 
see e.g. Ref. [1]. In the present case, it could be argued 
to fail because we are collapsing data at temperatures 
which are both above and below Tc. Indeed, below Tc 
scaling in the /5-regime is not expected anymore because 
the height of the plateau, fq in Eq. now becomes a 
temperature dependent quantity, with the following pre- 



IV. DYNAMIC HETEROGENEITY 

Having established the ability of MC simulations to 
efficiently reproduce the averaged slow dynamical be- 
haviour observed in ND simulations, we now turn to the 
study of the dynamic fluctuations around the average dy- 
namical behaviour, i.e. to dynamic heterogeneity. 

Dynamic fluctuations can be studied through a four- 
point susceptibility, Xi{t)i which quantifles the strength 
of the spontaneous fluctuations around the average dy- 
namics by their variance. 



X4(t) = iV„[(/2(q,t))-F2(g,t)] 



(13) 



where 



fs{ci,t) = — ^cos(q- [r,(t) ~r,(0)]), (14) 



represents the real part of the instantaneous value of the 
self-intermediate scattering function, so that Fa{q,t) — 
{fs{q_,t)). As shown by Eq. (fT^ . Xii^) will be large 
if run-to-run fluctuations of the self-intermediate scat- 
tering functions averaged in large but finite volume, are 
large. This is the case when the local dynamics becomes 
spati ally correlated, as already discussed in several pa- 
pers m, m, H, m, h 113, m . 

What X'i{t) captures is information on the spatial 
structure of the spontaneous fluctuations of the dynamics 
around their average. We define /i(q, t) = cos(q- [ri{t) — 
Ti (0)] ) , the contribution of particle i to the instantaneous 
value of Fg {q, t) , and 



SMq,t) ^ Mq,t) - F,iq,t), 



(15) 



10 




FIG. 9: Snapshots of dynamic heterogeneity at T = 6100, 
3580 and 3000 K (from top to bottom). The snapshot presents 
particles which, in a particular run at a particular tempera- 
ture have been slower than the average, and have therefore 
large, positive values of Sfi{q,t « r^) defined in Eq. (|15|l . 
Light colour is used for Si, dark for O. Slow particles tend 
to cluster in space on increasingly larger lengthscales when T 
decreases. 

its fluctuating part. Then Xii^) can be rewritten in the 
suggestive form, 

Xi{t)=p J rfr^|^(S/,(q,t)5/,(q,t),5(r-[r,(0)-r,(0)]) 

(16) 

where subtleties related to the exchange between thermo- 
dynamic limit and thermal average are discussed below. 



Therefore X4(i) is the volume integral of the spatial cor- 
relator between local fluctuations of the dynamical be- 
haviour of the particles. It gets larger when the spatial 
range of these correlations increases. 

To get a feeling of how these fluctuations look like in 
real space, we present snapshots at different tempera- 
tures in Fig. [21 To build these snapshots we show, for 
a given run at a given temperature, those particles for 
which the fluctuating quantity Sfi{q,t « Tq,), is posi- 
tive and larger than a given threshold, which wc choose 
close to 1 /2 for graphical convenience (this leads to about 
1/3 of the particles being shown, and clearer snapshots). 
The shown particles are therefore slower than the aver- 
age for this particular run. The evolution of the snap- 
shots between 6100 K and 3000 K clearly reveals the ten- 
dency for slow particles to cluster in space, revealing the 
growth of the lengthscale of kinetic heterogeneities. We 
should note, however, that the clusters shown here are 
not macroscopic objects even at the lowest temperature 
studied. Moreover, similar snapshots in Lennard- Jones 
systems reveal more clearly the tendency we seek to il- 
lustrate d^. We interpret this as a further qualitative 
indication that dynamic heterogeneity is less pronounced 
in this Arrheniusly relaxing material than in more fragile 
Lennard-Jones systems. 

We turn to more quantitative measures of dynamic het- 
erogeneity and show the time dependence of the dynamic 
susceptibility X4(0 obtained from our MC simulations 
for various temperatures in Fig. [TOl Similar data are 
obtained for Si and O, and we on ly p resent the former. 
As predicted theoretically in Ref. [43] we find that Xiit) 
presents a complex time evolution, closely related to the 
time evolution of the self-intermediate scattering func- 
tion. Overall, Xiit) is small at both small and large times 
when dynamic fluctuations are small. There is there- 
fore a clear maximum observed for times comparable to 
Tq,, where fluctuations are most pronounced. The posi- 
tion of the maximum then shifts to larger times when 
temperature is decreased, tracking the alpha-relaxation 
timescale. The most important physical information re- 
vealed by these curves is the fact that the amplitude of 
the peak grows when the temperature decreases. This is 
direct evidence, recall Eq. ([T6|), that spatial correlations 
grow when the glass transition is approached. 

The two-step decay of the self-intermediate scattering 
function translates into a two-power law regime for Xi (t) 
approaching its maximum. We have fitted these power 
laws, X4:{t) ^ followed by Xi{t) ^ t'' with the expo- 
nents a ~ 0.3 and b ~ 0.92 in Fig. [TUj We have inten- 
tionally used the notation a and b for these exponents 
which arc predicted, within mode-coupling theory, to be 
equal to the standard exponents also describing the time 
dependence of intermediate scattering functions [3, [13 ■ 
, Our findings are in reasonable agreement with values for 
a and b discussed above, although the 6-value is about 
50% too large. Moreover, a two-power law regime is only 
observed for T < Tc, where MOT does not apply any- 
more. We note that the 5-value is predicted to be 6 = 2 
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FIG. 10: Top: Four-point susceptibility, Eq. (fT3t . for the 
same temperatures as in Fig. O decreasing from left to right. 
The low temperature data at T = 2750 K are fitted with two 
power laws shown as dashed lines with exponents 0.3 and 0.92 
at short and large times, respectively. The envelope of the 
maxima is fit with an exponent 0.285. Bottom: temperature 
evolution of the maxima in various dynamic susceptibihties. 



from the perspective of modelling strong glass-formers 
using kinetically constrained models with Arrhenius be- 
haviour 49|; this prediction is clearly incorrect for BKS 
silica ji.[47f. 

We then focus on the amplitude of the dynamic suscep- 
tibility at its maximum and follow its temperature evo- 
lution in Fig. [ini As suggested by the snapshots shown 
in Fig. [9l we confirm that spatial correlations increase 
when T decreases, as X4 gets larger at low temperatures. 
The temperature evolution of the peak was discussed in 
Ref. [J. Both MCT and kinetically constrained models 
strongly overestimate the temperature evolution of X4 at 
its peak value, as emphasized already in Ref. Finally, 
we note that the typical values observed for the peak 
of X4 at low temperatures are significantly smaller than 
those observed for more fragile Lennard- Jones systems, 
suggesting once more that dynamic heterogeneity is less 
pronounced in strong glass-forming materials. 

This comparison is also useful to discuss the possibility 
of finite size effects on the present X4 data. If computed 
in a simulation box which is too small, the dynamic sus- 



ceptibility takes values that are too small [50| . Our data 
indicate that no saturation of the maximum value of X4 (t) 
is reached, and the values we find of smaller than the ones 
found in a Lennard- Jones system with a comparable sys- 
tem size and for which a detailed search for possible finite 
size effects was performed @. We believe therefore that 
our results are not affected by finite size effects. 

We then compare these results to the ones obtained 
using Newtonian dynamics in the same system. In that 
case, care must be taken of the order at which the ther- 
modynamic limit and the thermal average are taken in 
Eq. (Uni). Indeed when ND is used, the dynamics strictly 
conserves the energy during the simulation and thermal 
averages are then performed in the microcanonical en- 
semble, and xl' is measured. To measure xT in the 
canonical ensemble for ND, an additional contribution 
must be added, which takes into account the amount of 
spontaneous fluctuations which are due to energy fluctu- 
ations jsij , 



rp2 

xfit) = - 

cv 



dF,{q,t) 
dT 



—XT{t), (17) 
Cv 



where cy is the constant volume specific heat expressed 
in ks units. The results for xj and xf obtained from 
ND, and the difference term in Eq. p7|) . are all presented 
in Fig. [ini We find that the MC results for X4 lie closer 
to the microcanonical results obtained from ND, while 
the canonical fluctuations are significantly larger, due to 
the large contribution of the right hand side in Eq. (|17|) . 
This is at first sight contrary to the intuition that MC 
simulations are thermostatted and should be a fair rep- 
resentation of canonical averages in ND. But this is not 
what happens. As discussed in Refs. [1, 0|, a major role 
is played by conservation laws for energy and density 
when dynamic fluctuations are measured. In the case of 
energy conservation the mechanism can be physically un- 
derstood as follows. For a rearrangement to take place in 
the liquid, the system has to locally cross an energy bar- 
rier. If dynamics conserves the energy, particles involved 
in the rearrangement must borrow energy to the neigh- 
boring particles. This 'cooperativity' might be unneces- 
sary if energy can be locally supplied to the particles by 
an external heat bath, as in MC simulations. Conserva- 
tion laws, therefore, might induce dynamic correlations 
between particles and dynamic fluctuations can be dif- 
ferent when changing from Newtonian, energy conserv- 
ing dynamics to a stochastic, thermostatted dynamics. 
With hindsight, this is not such a surprising result. The 
specific heat, after all, also behaves differently in differ- 
ent statistical ensembles. The ensemble dependence and 
dependence upon the microscopic dynamics of dynamic 
susceptibilities in supercooled liquids are the main sub- 
jects of two recent papers [3, 0]. Our results for silica 
quantitatively agree with the theoretical analysis they 
contain, and with the corresponding numerical results 
obtained in Lennard-Jones systems. 

There is an experimentally extremely relevant conse- 
quence of these findings [5l|, [s^ . As shown in Fig. [TOl 



12 



the difference between the microcanonical and canoni- 
cal values of the dynamic fluctuations in ND represents 
in fact the major contribution to xj meaning that the 
term xf can be safely neglected in Eq. pT]) . Since the 
right hand side of pT]) is more easily accessible in an ex- 
periment than Xi itself, Eq. p7|) opens the possibility of 
an experimental estimate of the four-point susceptibility. 
This finding, and its experimental application to super- 
cooled glycerol and hard sphere colloids, constitute the 
central result of Ref. [5l| , while more data are presented 
in Ref. 

V. CONCLUSION 

We have implemented a standard Monte Carlo algo- 
rithm to study the slow dynamics of the well-known 
BKS model for silica in the temperature range from 
6100 K to 2750 K. Our results clearly establish that 
Monte Carlo simulations can be used to study the dy- 
namics of silica because quantitative agreement is found 
with results from Newtonian dynamics for the same po- 
tential, apart at very short times where thermal vibra- 
tions are efficiently suppressed by the Monte Carlo al- 
gorithm. The agreement between the two dynamics is 
by no means trivial and constitutes an important result 
of the present study. This suggests that Monte Carlo 
simulations constitute a useful and efficient tool to study 
also the nonequilibrium aging dynamics of glass- forming 
liquids, a line of research initiated in Ref. [53|. 

Since dynamical correlations arc not affected by short- 
time vibrations, we have been able to revisit the mode- 
coupling analysis initially performed in Ref. p^ . We 
find that mode-coupling theory accounts for the quali- 
tative features of the data quite well, but the detailed, 
quantitative predictions made by the theory were shown 
to fail: correlation functions close to the plateau do not 
follow the behaviour predicted for the MCT /3-correlator, 
time-temperature superposition only holds at very large 
times but fails at smaller times because the plateau in 
correlation function is strongly temperature dependent, a 
dependence which does not follow the singular behaviour 
predicted by MCT. Moreover, the temperature regime 
where the theory can supposedly be applied is found to 
be at most 1 decade when only the temperature evolu- 
tion of relaxation timescales is considered. Furthermore, 
we have argued that the motivation to analyze silica data 
in terms of MCT, a "fragile to strong" crossover, can in 



fact be more simply accounted for in terms of a crossover 
between two distinct Arrhenius regimes occurring close 
to T « 4000 K. Overall these results suggest a negative 
answer to the question: is there any convincing evidence 
of an avoided mode-coupling singularity in silica? 

We have finally analyzed dynamic heterogeneity in sil- 
ica. Wc find that the dynamics is indeed spatially het- 
erogeneous, and spatial correlation of the local dynami- 
cal behaviour was shown to increase when temperature 
decreases. Wc also found that all indicators of dynam- 
ically heterogeneous dynamics such as decoupling and 
four-point dynamic susceptibilities, suggest that the ef- 
fects are less pronounced in silica than in more fragile 
glass-forming materials, but do not seem qualitatively 
different. 

The most natural interpretation is that strong and 
fragile materials in fact belong to the same class of ma- 
terials, where the effects of dynamic heterogeneity could 
become less pronounced, but definitely non-zero, for ma- 
terials with lower fragility. This suggests that it could 
be incorrect to assume that strong materials belong to 
a different universality class from fragile ones, as studies 
of kinctically constrained models with different fragili- 
ties would suggest [i^ [131 , and they should rather seat 
at the end of the spectrum of fragile systems. It seems 
however similarly incorrect to consider that strong ma- 
terials are "trivial" because an Arrhenius behaviour can 
be explained from simple thermal activation over a fixed 
energy barrier corresponding to a local, non-cooperative 
event. Our results show that this is not a correct rep- 
resentation of the physics of strong glass-formers either. 
Convincingly incorporating fragility into current theories 
of the glass transition while simultaneously giving it a mi- 
croscopic interpretation remains therefore an important 
challenge. 
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